This notebook provides functions to calculate the transport of water through packages and calculates the partitioning of water between water vapor and solid materials. The calculations can be used to simulate the time dependent concentration of water in various portions of a pharamaceutical packaging system.
In [52]:
import numpy as np
from datetime import datetime, timedelta
import pandas as pd
import argparse as ap
from matplotlib import pyplot as plt
from scipy.optimize import minimize_scalar
from scipy import interpolate
from pdb import set_trace
import json
%matplotlib inline
In [53]:
store = pd.HDFStore("simulation_constants.hdf", mode="r")
GAB_constants = store["GAB_constants"]
package_properties = store["package_properties"]
store.close()
In [54]:
package = {
"start_date": "2016-04-01",
"end_date": "2017-04-01",
"rh": "40",
"temperature": "30",
"layers": [{
"container": {
"ID": "BOTTLE_45CC",
"seal": "UNACTIVATED"
},
"desiccant": {
"type": "SILICA",
"mass": "10",
"water_content": "20"
},
"product": {
"units": "100",
"unit_mass": "20.",
"bulk_density": "1.",
"components": [
{"name": "GENERIC_API",
"frac": "0.1",},
{"name": "MANNITOL",
"frac": "0.3"},
{"name": "LACTOSE_REGULAR",
"frac": "0.6"}
]
}
}],
"events": [{
"date": "2016-04-01",
"layer": "0",
"ID": "replace_desiccant",
"desiccantWater": "0.2",
"frequency": "WEEKLY"
},
{
"date": "2016-04-03",
"layer": "0",
"ID": "remove_product",
"units": 1,
"frequency": "DAILY"
}
]
}
jsystem = json.dumps(system)
In [83]:
class Container(object):
"""
Define the type of the container and necessary information, such as
volume, area, ... used for calculating transport through the container material.
Class variables:
ID : string - unique identification string to lookup parameters
name : string - name of the container (e.g. LDPE)
volume : float - total volume of container (cm^3)
area : float - surface area of package (cm^2)
thickness : float - thickness of container used for flux calculation (mil)
Only used for LDPE.
seal : string - key work identifier describing the permeability of the seal
transport_coef : dictionary {"slope": float, "intercept": float} - used for
calculating the permeability of the package.
temperature : float - container environment temperature (K)
permeability : permeability of the container at the set temperature (mg / day / delta RH)
Class methods:
set_properties : Set the properties of the container.
recalc_properties : Calculate and set all the container properties based on the current
object variable values.
set_temperature : Set the temperature of the container.
set_seal : Set the type of seal for the container.
calc_area : Calculate the area of the container
calc_permeability : Calculate the permeability
"""
def __init__(self, ID, **kwargs):
"""
Set the class attributes.
"""
self.set_properties(ID, **kwargs)
def set_properties(self, ID, temperature=298.15, **kwargs):
"""
Set all the properties of the package.
Parameters:
ID: string - Unique identification of the container
temperature: float - temperature of environment (K)
optional kwargs:
seal : string - value for the seal type
volume : float - volume of the container
thickness : float - thickness for LDPE
area : float - area of LDPE
temperaure : float - temperature of container environment
"""
store = pd.HDFStore("simulation_constants.hdf", mode="r")
package_properties = store["package_properties"]
store.close()
#set_trace()
if ID in package_properties.index.values:
self.ID = ID
self.name = package_properties.loc[ID]["name"]
self.transport_coef = package_properties.loc[ID][["slope", "intercept"]].to_dict()
if ID != "LDPE":
self.volume = package_properties.loc[ID]["volume"]
else:
raise ValueError("No container with identification %s is known."%ID)
if ID!="LDPE":
if "seal" in kwargs:
self.seal = kwargs["seal"]
else:
self.seal = "ACTIVATED"
elif ID=="LDPE":
if "volume" not in kwargs or thickness not in kwargs:
raise ValueError("Must supply volume and thickness for LDPE.")
else:
self.volume = float(volume)
self.thickness = float(thickness)
if "area" in kwargs:
self.area = float(kwargs["area"])
else:
self.area = calc_area()
self.temperature = float(temperature)
self.permeability = self.calc_permeability()
def recalc_properties(self):
"""
Recalculate and set all the container properties based on the
new object variable values.
"""
self.permeability = self.calc_permeability()
def set_temperature(self, temperature):
"""
Set the environment temperature for the container and recalculate
container parameters.
Parameters
temperature: float - temperature of the environment (K)
"""
self.temperature = float(temperature)
self.recalc_properties()
def set_seal(self, seal):
"""
Set the seal conditions and recalculate the container parameters.
Parameters
seal : string - type of seals
["ACTIVATED" | "UNACTIVATED" | "BROACHED"]
"""
if seal in ["ACTIVATED", "UNACTIVATED", "BROACHED"]:
self.seal = seal
else:
raise ValueError("Unknown seal type")
self.recalc_properties()
def calc_area(self):
"""
Calculate the package area assuming the volume is a sphere.
return : float - Area of package in cm^2.
"""
area = 4 * np.pi * (3 * self.volume / (4 * np.pi))**(2./3.)
return area
def calc_permeability(self):
"""
Calculate the permeability of the package.
return : float - permeabilty mg / day / delta RH
"""
if self.ID == "LDPE":
permeability = np.exp(self.transport_coef["intercept"] \
+ self.transport_coef["slope"] / self.temperature) \
* self.area * 4. / (10. * self.thickness)
else:
if self.seal == "UNACTIVATED":
f1 = 1.20923913
f2 = 0.
elif self.seal == "BROACHED":
f1 = 1.
f2 = 0.02072
else:
f1 = 0.
f2 = 0.
permeability = np.exp(self.transport_coef["intercept"] \
+ self.transport_coef["slope"] / self.temperature) * f1 + f2
return permeability
@classmethod
def make_container(cls, config):
"""
Make the container object using the dictionary for parameters.
config : dictionary - Dictionary of parameters that define the
container.
"""
ID = config["ID"]
del config["ID"]
container = cls(ID, **config)
return container
In [56]:
class Desiccant(object):
"""
Define the desiccant inside the container (e.g. type, amount, water...etc).
Class Attributes
ID : string - unique identification number to lookup parameters
name : string - Desiccant material.
mass : float - mass of desiccant (g)
water_content : float - mass fraction of water contained in the desiccant
(mass water (mg) / mass dry desiccant (g))
GAB_parameters : dictionary - {Wm, C, K} GAB constants from lookup file.
density: float - density of dry desiccant (g/cm^3)
water : float - total mass of water (g)
Class Methods
set_properties : Set the properties of the desiccant
refresh : Refresh the desiccant to a new water content
equilibrate : Equilibrate the deisccant to a provided RH and set the water content.
calc_water_content: Calculate the water content from GAB and passed water activity
"""
def __init__(self, ID, mass, **kwargs):
self.set_properties(ID, mass, **kwargs)
def set_properties(self, ID, mass, density=1.0, **kwargs):
"""
Set the properties of the desiccant.
Parameters
ID : string - Unique identification of the desiccant for the lookup.
mass : float - Mass of the desiccant.
optional kwargs
water_content : float - mass fraction of water in desiccant
(mass water (mg) / mass dry desiccant (g))
density : float - density of dry desiccant (g)
"""
store = pd.HDFStore("simulation_constants.hdf", mode="r")
GAB_constants = store["GAB_constants"]
store.close()
if ID in GAB_constants.index.values:
self.ID = ID
self.name = GAB_constants.loc[material]["name"]
self.GAB_parameters = GAB_constants.loc[ID][["C","K","Wm"]].to_dict()
else:
raise ValueError("Desiccant type %s is not defined" %ID)
self.mass = float(mass)
if water_content in kwargs:
self.water_content = float(kwargs["water_content"])
elif initial_water_activity in kwargs:
self.water_content = \
self.calc_water_content(kwargs["initial_water_activity"]) * 1.e3
else:
self.water_content = 20.
self.density = float(density)
self.water = self.water_content * self.mass * 1.e-3
def refresh(self, water_content=20., initial_activity=None):
"""
Refresh the desiccant (e.g. replace with equivalent desiccant
mass with lower water content). Specify either new water
content (water_content), or initial water activity of the
desiccant (initial_activity).
Parameters
water_content: float - Water content of the fresh desiccant
(mg water / g dry desiccant)
initial_activity: float - Water activity of the fresh desiccant (unitless)
"""
if initial_activity is None:
self.water_content = float(water_content)
else:
self.water_content = \
self.calc_water_content(float(initial_activity))
self.water = self.water_content * self.mass * 1.e-3
def equilibrate(self, aw):
"""
Equilibrate the desiccant to a new water activity.
Parameters
aw: float - Water activity to equilibrate desiccant.
"""
self.water_content = self.calc_water_content(aw)
self.water = self.water_content * self.mass * 1.e-3
def calc_water_content(self, aw):
"""
Calculate the water content from the GAB parameters at the provided
water activity.
Parameters
aw : float - Water activity
return: water_content (mg water / g desiccant)
"""
aw = float(aw)
Wm = self.GAB_parameters["Wm"]
C = self.GAB_parameters["C"]
K = self.GAB_parameters["K"]
water_content = (Wm*C*K*aw) / ((1-K*aw) * (1-K*aw+C*K*aw)) * 1.e3
return water_content
In [57]:
class Product(object):
"""
Define the product contents inside the container (e.g. units,
unit mass, unit composition...etc.). Provides calculation of the
average GAB parameters over the individual components.
Class Attributes:
name: string - name of the product
units: integer - number of unit deses of product
unit_mass: float - mass of a single unit (mg/unit)
unit_volume: float - volume of a single unit (cm^3/unit)
bulk_density: float - density of the bulk substance (g/cm^3)
components: list of dictionaries [{"ID": string, "frac": float,
"GAB_parameters": dict}]
water_content: float - mass fraction of water in product (unitless)
(mass water (g) / mass dry product (g))
avg_GAB: dictionary of average GAB parameters for the product.
mass: float - total mass of the product (g)
water : float - total mass of water (g)
Class Methods
set_properties: set the properties of the product
(number units, components, ...etc)
remove_product: Remove some units/mass of the product and recalculate units/mass.
equilibrate : Equilibrate the product to a new RH and set the new water content.
calc_water_content: Calculate the water content of the product for a given aw.
calc_avg_GAB: Calculate the average GAB parameter for the bulk composition.
calc_mass: Calculate the total product mass
"""
def __init__(self, units, unit_mass, unit_volume, components, **kwargs):
self.set_properties(units, unit_mass, unit_volume, components, **kwargs)
def set_properties(self, units, unit_mass, unit_volume, components,
bulk_density=1.0, **kwargs):
"""
Set the properties of the product.
Parameters:
name: string - Name of product
unit_mass: float - mass of single unit (mg/unit)
unit_volume: float - volume of single unit (cm^3/unit)
units: integer - number of unit doses
components: list of dictionaries [{"ID"}: string, "frac": float,
"GAB_parameters": dict}]
Optional Parameters
bulk_density: float - bulk density of the product (default 1 g/cm^3)
water_content: float - water content of the product (unitless)
(mass water / mass dry product)
initial_water_activity: float - initial water activity to calculate
the water content from (unitless)
"""
store = pd.HDFStore("simulation_constants.hdf", mode="r")
GAB_constants = store["GAB"]
store.close()
if "name" in kwargs:
self.name = kwargs["name"]
self.units = units
self.unit_mass = unit_mass
self.unit_volume = unit_volume
self.mass = self.units * self.unit_mass * 1.e-3
for c in components:
if c["ID"] in GAB_constants.index.values:
c["GAB_parameters"] = GAB_constants.loc[ID][["C","K","Wm"]].to_dict()
else:
raise ValueError("Product type %s is not defined" %c["ID"])
self.avg_GAB = self.calc_avg_GAB()
self.bulk_density = float(bulk_density)
if "water_content" in kwargs:
self.water_content = float(kwargs["water_content"])
elif "initial_water_activity" in kwargs:
self.water_content = \
self.calc_water_content(kwargs["initial_water_activity"])
else:
self.water_content = self.calc_water_content(0.4)
self.water = self.water_content * self.mass
def remove_product(self, units=None, mass=None):
"""
Reduce the mass of the desiccant. Specify either units or
mass.
Parameters
units: float - Number of units to remove.
mass: float - Mass to remove (g)
"""
if units is not None:
self.units = self.units - units
self.mass = self.units * self.unit_mass * 1.e-3
elif units is None and mass is not None:
self.mass = self.mass - mass
self.units = self.mass * 1.e3 / self.unit_mass
else:
raise ValueError("Must supply either units or mass.")
self.water = self.water_content * self.mass
def equilibrate(self, aw):
"""
Equilibrate the product to a new water activity.
Parameters
aw: float - Water activity to equilibrate product.
"""
self.water_content = self.calc_water_content(aw)
self.water = self.water_content * self.mass
def calc_avg_GAB(self):
"""
Calculate the average GAB parameters for the composition of components.
See "GAB Generalized Equation for Sorption Phenomena", Food and Bioprocess Technology
March 2008 Vol 1, Issue 1, pp 82--90
aw/w = a + b*aw + c*aw^2
a = 1/(Wm*C*K)
b = (C-2)/(Wm*C)
c = K(1-C)/(Wm*C)
"""
def GAB(Wm, C, K, aw):
return (Wm*C*K*aw) / ((1-K*aw) * (1-K*aw+C*K*aw))
aw = np.arange(0.1, 0.8, 0.05)
w = np.array([np.sum([c["frac"] * GAB(c["Wm"], c["C"], c["K"], aw_i)
for c in self.components]) for aw_i in aw])
y = aw / w
[c, b, a] = np.polyfit(aw, y, 2)
K = (-b + np.sqrt(b**2 - 4*a*c))/(2*a)
C = b / (a*K) + 2
Wm = 1 / (b + 2*K*a)
return {"K": K, "C":C, "Wm":Wm}
def calc_water_content(self, aw):
"""
Calculate the water content of the product mixture at activity, aw, based on the
Guggenheim, Anderson and de Boer (GAB) three-parameter isotherm model.
See "GAB Generalized Equation for Sorption Phenomena", Food and Bioprocess Technology
March 2008 Vol 1, Issue 1, pp 82--90
w = (Wm*C*K*aw) / ((1-K*aw)*(1-K*aw+C*K*aw))
Parameters:
aw: float - water activity
return float - mass fraction of water (unitless) (mass of water / mass of dry substance)
"""
Wm = self.avg_GAB["Wm"]
C = self.avg_GAB["C"]
K = self.avg_GAB["K"]
water_content = (Wm*C*K*aw) / ((1-K*aw) * (1-K*aw+C*K*aw))
return water_content
In [58]:
class Layer(object):
"""
Define a layer of the package, which consists of a container with desiccant and
product.
Class Attributes:
container: object - instance of Container class
desiccant: object - instance of Desiccant class
product: object - instance of Product class
temperature: float - temperature (K)
rh: float - relative humidity (unitless)
head_space: float - current head space in the layer (cm^3)
isotherm: function - isotherm of total layer
function signature:
isotherm(water_mass) -> aw
water: float - total mass of water in the system (vapor, product, desiccant) (g)
Class methods
set_properties : Calculate and set the layer properties based on input arguments.
remove_product : Remove product from the product child and recalculate properties.
replace_deisscant : Replace the desiccant with fresh desiccant and recalculate properties.
calc_isotherm : Calculate the isotherm of the layer
equilibrate : Equilibrate the layer as a closed isothermal system.
calc_head_space : Calculate the head space of the layer
calc_water_mass : Calculate the mass of water in the layer.
"""
def __init__(self, container, desiccant, product, temperature, rh, **kwargs):
self.set_properties(container, desiccant, product, temperature, rh, **kwargs)
def set_properties(self, container, desiccant, product, temperature, rh, **kwargs):
"""
Set the properties of the layer
"""
self.container = container
self.desiccant = desiccant
self.product = product
self.temperature = temperature
self.rh = rh
self.head_space = self.calc_head_space()
self.water = self.calc_water_mass()
self.isotherm = self.calc_isotherm()
self.equilibrate()
def remove_product(self, units=None, mass=None):
"""
Remove some product from the layer. Recalculates
the amount of product, the new head space and the
new isotherm of the layer.
Parameters
units: float - Number of units of product to remove (dimensionless)
"""
# First calculate the new product units
self.product.remove_product(units=units, mass=mass)
# Next calculate the new head space.
self.head_space = self.calc_head_space()
# Now calculate the new total wate mass
self.water = self.calc_water_mass()
# Finally, calculate the new isotherm
self.isotherm = self.calc_isotherm()
# Note that the RH does not change.
def replace_desiccant(self, water_content=20., initial_activity=None):
"""
Replace the desiccant with fresh desiccant. Recalculates
the new isotherm of the layer. If initial_activity is provided,
uses initial_activity, else uses water_content (default 20 mg/g).
Parameters
water_content: float - The water content of the fresh deisccant
(mg water / g dry desiccant)
default 20 mg / g
initial_activity: float - The initial water activity of the fresh
desiccant (unitless)
"""
# First set the new desiccant water content
self.desiccant.refresh(water_content, initial_activity)
# Next caluclate the new mass of water in the layer.
self.water = self.calc_water_mass()
# Finally, calculate the new RH (Note that the isotherm does not change.).
self.rh = self.equilibrate()
def calc_head_space(self):
"""
Calculate the head space of the layer.
return float - head space of layer (cm^3)
"""
volume_total = self.container.volume
volume_desiccant = self.desiccant.mass / self.desiccant.density
volume_product = self.product.mass / self.product.bulk_density
head_space = volume_total - volume_desiccant - volume_product
return head_space
def calc_isotherm(self):
"""
Calculate the effective isotherm of the layer
(product + desiccant + head space).
return: function - isotherm of the layer
"""
aw = np.linspace(0.1, 0.8, 100)
T = self.T
dens_sat = 5.079e-3 + 2.5547e-4*T + 1.6124e-5*T**2 + 3.6608e-9*T**3 \
+ 3.9911e-9*T**4
water = np.array([dens_sat * self.head_space() * 1.e-3 * aw_i \
+ self.desiccant.calc_water_content(aw_i) * self.desiccant.mass * 1.e-3 \
+ self.product.calc_water_content(aw_i) * self.product.mass \
for aw_i in aw])
isotherm = interpolate.interp1d(water, aw)
return isotherm
def equilibrate(self):
"""
Equilibrate the system as a closed system for mass. That is
calculate the equilibrium distribution of water keeping the
total water mass constant.
return float - water activity (unitless)
"""
# First calculate the equilibrium water activity (RH).
aw = self.isotherm(self.water_mass)
# Now calculate the new desiccant and product water content.
self.desiccant.equilibrate(aw)
self.product.equilibrate(aw)
def calc_water_mass(self):
"""
Calculate the total mass of water in the layer (head space,
product and desiccant)
return: float - mass of water (g)
"""
T = self.T
m_v = 5.079e-3 + 2.5547e-4*T + 1.6124e-5*T**2 + 3.6608e-9*T**3 \
+ 3.9911e-9*T**4 * self.rh * self.head_space * 1.e-3
mass = m_v + self.desiccant.water + self.product.water
return mass
def open_layer(self, ambient_rh):
"""
Open the container of the layer, exposing to ambient conditions.
Then recalculate the new mass of water and equilibrate.
Parameters
ambient_rh : float - Ambient relative humidity the open container is
exposed. (unitless)
"""
self.rh = ambient_rh
self.water = self.calc_water_mass()
self.equilibrate()
def flux(self, ambient_rh, time_step):
"""
Calculate the mass of water that permeates from the container
exterior into the layer. Then calculate the new mass of water
and equilibrate.
Parameters
ambient_rh : float - The relative humidity on the outside of the
container (unitless).
time_step : float - The duration of time used to calculate the total
mass that permeates through the container (days).
"""
delta_rh = (ambient_rh - self.rh)
delta_water = self.container.permeability * delta_rh * time_step
self.water += delta_water
self.equilbrate()
In [59]:
class Events(object):
"""
Define and events object itterator that returns events in
chronological order. This is non-trivial as the events
have a defined date and can have a reoccurring frequency.
Class attributes
events : list of dictionaries - {date, ID, frequency, [mass, water]}
current_date : datetime object - Defines the most current
event date that the generator has provided.
last_date : datetime object - Defines the last date to generate
events.
"""
def __init__(self, events, date_format="%Y-%m-%d",
last_possible_date=datetime(3000, 1, 1)):
"""
Convert the date string in each event to a datetime
object, and add a timedelta key-value pair to the
dictionary that relates to the reoccurrence frequency.
Parameters
events : list of dictionaries
date_format : string - Defines the date format of the
string contained in the event.
last_possible_date : datetime object - latest possible date
for retrieving event.
"""
last_date = datetime.strptime(events[0]["date"], date_format)
current_date = datetime.strptime(events[0]["date"], date_format)
for e in events:
# Convert the string date to a datetime object.
e["date"] = datetime.strptime(e["date"], date_format)
# Set the current date to the earliest event date.
if e["date"] < current_date:
current_date = e["date"]
# Set the last date to either the last explicit date,
# if no event is reoccuring, or to an arbitrary VERY distant date.
if e["date"] > last_date:
last_date = e["date"]
# Add to the dictionary the time delta corresponding
# to the frequency.
if e["frequency"] == "ONCE":
# For one time events, just use a VERY large interval
e["time_delta"] = datetime.timedelta(99999)
elif e["frequency"] == "DAILY":
e["time_delta"] = datetime.timedelta(1)
last_date = last_possible_date
elif e["frequency"] == "WEEKLY":
e["time_delta"] = datetime.timedelta(7)
last_date = last_possible_date
elif e["frequency"] == "MONTHLY":
e["time_delta"] = datetime.timedelta(30)
last_date = last_possible_date
else:
raise ValueError("Event frequency %s is not known" \
%e["frequency"])
self.current_date = current_date
self.last_date = last_date
self.events = events
def __iter__(self):
return self
def __next__(self):
"""
Returns the next chronological events.
returns : list - List of events that are in next
chronological order.
"""
current_events = []
while len(current_events) == 0:
for e in self.events:
# Include any event that is on a multiple of it's
# reoccurence interval (including 0).
# Note that mod(0,#)=0 and mod(#,0)=0.
if np.mod((self.current_date - e["date"]).days, \
e["time_delta"].days) == 0:
current_events.append(e)
self.current_date += timedelta(1)
return current_events
In [65]:
class Package(object):
"""
Define the entire package system, consisting of an ambient
conditions, multiple layers, and events.
Class attributes
layers : list of Layer objects - list of the layers in the package
Note that the outer most layer should be the first item in
the list and the inner most layer the last item in the list.
temperature : float - ambient temperature (C)
rh : float - ambient relative humidity (unitless)
start_date : datetime - start date for the simulation
end_date : datetime - end date for the simulation
events : Event object
results : list
Class methods
simulate_interval : Simulate the package dynamics for a time interval
that the package configuration does not change.
execute_event : Execute an event defined by the user.
simulate : The main function of the class to perform the simulation
of the package over the entire input time with possible
intermediate events.
"""
def __init__(self, layers, temperature, rh, start_date, end_date,
events=None):
self.layers = layers
self.temperature = temperature
self.rh = rh
self.start_date = datetime.strptime(start_date, "%Y-%m-%d")
self.end_date = datetime.strptime(end_date, "%Y-%m-%d")
self.event = events
# Initialize the results with the beginning configuration.
self.results = [0, [l.rh for l in layers]]
def simulate_interval(self, duration):
"""
Simulate the package dynamics from one event
to the next event.
Parameters
duration : float - total time to perform the calculation
(elapsed time between events.)
"""
# Set the starting time to the previous starting time.
t_o = results[-1][0]
t = 0
# Begin loop dynamics
while t < duration:
# Set the time step
dt = 1.e-10
t += dt
# Calculate the flux through the individual layers.
self.layers[0].flux(self.rh, dt)
for i in range(1, len(self.layers)):
self.layers[i].flux(self.layers[i-1], dt)
# Equilibrate the layers.
for layer in self.layers:
layer.equilibrate()
self.results.append([t_o+t, [l.rh for l in self.layers]])
def simulate(self):
"""
Simulate the package dynamics. Creat a list of events from
the event object, perform the event then call the simulate_event
method to calculate the dynamics between events.
"""
if self.events is None:
# For no events, simply simulate the package for
# the requested duration.
duration = (self.end_date - self.start_date).days
self.simulate_interval(duration)
else:
# Get a chronological list of the events.
event_list = self.events.make_event_list()
# Set the current time point to the overall simulation
# start date.
current_date = self.start_date
# Cycle through the events.
for e in event_list:
if (e["date"] - current_date).days < 1:
# If the interval between the current date and
# the next event date is less than 1 day, then
# modify the package configuration before beginning
# the next simulation.
execute_event(e)
else:
# If the interval is greater than 1 day, then
# run a simulation for the interval.
duration = (e["date"] - current_date).days
self.simulate_interval(duration)
# Set the new current date to the last run event date.
current_date = e["date"]
def execute_event(self, event):
"""
Modify the package configuration according to the
event.
Parameters
event : dictionary - Dictionary that describes the event
{date, ID, [additional arguments]}
"""
# Get the command and layer number.
cmd = event["ID"]
layer_no = event["layer"]
# Remove the command and layer number from the dictionary.
del event["ID"]
del event["layer"]
# Execute the command on the layer and pass the remaining
# dictionary.
self.layers[layer_no].__getattribute__(cmd)(**event)
@classmethod
def make_package(cls, config):
"""
Create the package, layers, product, desiccant and container
objects from the config strgin.
Parameters
config : dictionary - Dictionary of all the parameters
required to build the package.
"""
for l in config["layers"]:
container = Container(l["container"]["ID"],
temperature=config["temperature"],
**l["container"])
@staticmethod
def convert_json(config):
"""
Convert the json string into a python dictionary. Attempt to
convert any stings that contain numberical values to floats
(or maybe integers).
"""
def convert(dictionary):
"""
This function is called recursively to try to
change the values of the items.
"""
new_dict = dict()
for k, v in dictionary.items():
if type(v) is dict:
new_dict[k] = convert(v)
elif type(v) is list:
new_dict[k] = list()
for itm in v:
if type(itm) is dict:
new_dict[k].append(convert(itm))
else:
try:
new_dict[k] = float(v)
except TypeError:
new_dict[k] = v
except ValueError:
new_dict[k] = v
return new_dict
converted_config = convert(config)
return converted_config
In [66]:
p = Package.make_package(package)
In [62]:
package = {
"start_date": "2016-04-01",
"end_date": "2017-04-01",
"rh": "40",
"temperature": "30",
"layers": [{
"container": {
"ID": "BOTTLE_45CC",
"seal": "UNACTIVATED"
},
"desiccant": {
"type": "SILICA",
"mass": "10",
"water_content": "20"
},
"product": {
"units": "100",
"unit_mass": "20.",
"bulk_density": "1.",
"components": [
{"name": "GENERIC_API",
"frac": "0.1",},
{"name": "MANNITOL",
"frac": "0.3"},
{"name": "LACTOSE_REGULAR",
"frac": "0.6"}
]
}
}],
"events": [{
"date": "2016-04-01",
"layer": "0",
"ID": "replace_desiccant",
"desiccantWater": "0.2",
"frequency": "WEEKLY"
},
{
"date": "2016-04-03",
"layer": "0",
"ID": "remove_product",
"units": 1,
"frequency": "DAILY"
}
]
}
jpackage = json.dumps(package)
In [133]:
class simulate(object):
def __init__(self, system):
"""
Extract the variables from the system JSON object
"""
self.json = system
self.system = json.loads(system)
self.store = pd.HDFStore("simulation_constants.hdf", mode="r")
def verify_system(self):
"""
Verify the system object make any format changes.
"""
def convert(dictionary):
new_dict = dict()
for k, v in dictionary.iteritems():
if type(v) is dict:
new_dict[k] = convert(v)
elif type(v) is list:
new_dict[k] = list()
for itm in v:
if type(itm) is dict:
new_dict[k].append(convert(itm))
else:
try:
new_dict[k] = float(v)
except TypeError:
new_dict[k] = v
except ValueError:
new_dict[k] = v
return new_dict
self.system = convert(self.system)
self.system["conditions"]["temperature"] += 273.
def set_bulk_density(self):
"""
Set a default density of the inner package, if it isn't supplied.
"""
if not self.system["inner"]["product"]["bulk_density"]:
self.system["inner"]["product"]["bulk_density"] = 1.3
def set_volume(self, loc="inner"):
"""
Set volume of the package, if it isn't set, by looking up in a reference table or
calculating it from the package contents.
"""
pkg = self.system[loc]
prod = pkg["product"]
if pkg["type"] == "LDPE":
if not pkg["volume"]:
pkg["volume"] = prod["units"] * ( prod["unit_mass"] / 1000. ) \
/ prod["bulk_density"]
if not pkg["area"]:
pkg["area"] = 4 * np.pi * (3 * pkg["volume"] / (4 * np.pi)) ** (2./3.)
else:
pkg["volume"] = self.store["volume"].loc[pkg["type"]]["volume"]
def set_desiccant_GAB(self, loc="inner"):
"""
Set the desiccant GAB parameters for the package.
"""
pkg = self.system[loc]
pkg["desiccant"]["GAB"] = \
self.store["GAB"].loc[pkg["desiccant"]["type"]][["Wm", "C", "K"]].to_dict()
def set_product_GAB(self, loc="inner"):
"""
Set the product average GAB parameters.
"""
pkg = self.system[loc]
components = pkg["product"]["components"]
frac = [np.float(c["frac"]) for c in components]
names = [c["name"] for c in components]
GABs = self.store["GAB"].loc[names][["Wm", "C", "K"]].to_dict("records")
pkg["product"]["GAB"] = GAB_avg(frac, GABs)
def set_transport_coef(self, loc="inner"):
"""
Set the water vaport transport coefficient for the package by looking it up
the parameters in a reference table.
"""
pkg = self.system[loc]
T = self.system["conditions"]["temperature"]
if pkg["type"] == "LDPE":
m, b = self.store["transport/LDPE"].iloc[0][["slope", "intercept"]]
A = pkg["area"]
l = pkg["thickness"] if pkg["thickness"] else 1.
transport_coef = transport_coef_calc(T, m, b, l=l, A=A, f1=1, f2=0)
else:
m, b = self.store["transport/package"].loc[pkg["type"]][["slope", "intercept"]]
if pkg["seal"] == "ACTIVATED":
f1, f2 = 1., 0.
elif pkg["seal"] == "UNACTIVATED":
f1, f2 = 1.21, 0.
elif pkg["seal"] == "BROACHED":
f1, f2 = 1., .0207
transport_coef = transport_coef_calc(T, m, b, l=1, A=1, f1=f1, f2=f2)
pkg["permeability"] = transport_coef
In [134]:
sim = simulate(jsystem)
sim.verify_system()
sim.set_volume(loc="inner")
sim.set_desiccant_GAB(loc="inner")
sim.set_product_GAB(loc="inner")
sim.set_transport_coef(loc="inner")
In [135]:
sim.system
Out[135]:
In [35]:
GAB_params = GAB_constants.loc[[component["name"] for component in
system["inner"]["product"]["components"]]][["Wm", "C", "K"]]
frac = np.array([component["frac"] for component in system["inner"]["product"]["components"]], dtype="f8")
In [36]:
# Define a few constants to use in the tests
TEMPERATURE = np.arange(280,340,5)
RH = np.arange(0.1, 0.9, 0.05)
WATER_ACTIVITY = np.arange(0.1,0.8,0.05)
In [1]:
def mass_density_sat(T):
"""
Mass of water in one cubic meter of saturated air at one bar and temperature T
parameters:
T: float - Temperature (K)
returns float - mass of water in one cubic meter saturated air (kg/m^3)
"""
T = T + 273.
return (5.079e-3 + 2.5547e-4*T + 1.6124e-5*T**2 + 3.6608e-9*T**3 + 3.9911e-9*T**4) / 1000.
In [38]:
plt.plot(TEMPERATURE, map(mass_density_sat, TEMPERATURE))
Out[38]:
In [39]:
def mass_water_vapor(T, rh, V):
"""
Calculate the total mass of water vapor in air at a
specified rh and a given size container.
parameters
T: float - Temperature (K)
rh: float - relative humidity (unitless)
V: float - volume of vapor (m^3)
returns float - mass of water in volume V (kg)
"""
return mass_density_sat(T) * rh * V
In [40]:
m = np.array([mass_water_vapor(323., rh_i, 1.e-4) for rh_i in RH])
plt.plot(RH, m)
Out[40]:
In [41]:
def GAB(aw, Wm, C, K):
"""
Calculate the water content of a substance based on the Guggenheim, Anderson and de Boer (GAB) three-parameter
isotherm model. See "GAB Generalized Equation for Sorption Phenomena", Food and Bioprocess Technology
March 2008 Vol 1, Issue 1, pp 82--90
w = (Wm*C*K*aw) / ((1-K*aw)*(1-K*aw+C*K*aw))
parameters:
aw: float or numpy array - water activity
Wm: float - GAB parameter
C: float - GAB parameter
K: float - GAB parameter
returns float - water content of substance (mass water substance (kg) / mass substance (kg))
"""
return (Wm * C * K * aw) / ((1 - K * aw) * (1 - K * aw + C * K * aw))
In [42]:
water = np.array([GAB(WATER_ACTIVITY, *GAB_params.loc[c].values) for c in GAB_params.index.values])
plt.plot(WATER_ACTIVITY, water.T)
plt.legend(GAB_params.index.values, loc="best")
Out[42]:
In [43]:
def mass_water_solid_mixture(aw, frac, params):
"""
Calculate the mass of water in a solid mixture at a specified water activity using
the superposition of GAB parameters.
parameters:
aw: float - water activity
frac: list of floats - list of mass fractions of the individual components [f_i] ; for i=1...N
params: list of dictionaries - list of GAB parameters [{wm, C, K}_i] ; for i=1...N
returns float - water content of substance (mass water substance (kg) / mass substance (kg))
"""
return np.sum([frac[i] * GAB(aw, p["Wm"], p["C"], p["K"]) for i, p in enumerate(params)])
In [44]:
avg = np.array([mass_water_solid_mixture(aw_i, frac, GAB_params.to_dict("records")) for aw_i in WATER_ACTIVITY])
plt.plot(WATER_ACTIVITY, avg)
Out[44]:
In [45]:
def GAB_regress(aw, w):
"""
Calculate the GAB parameters from water content - humidity measurements.
See "GAB Generalized Equation for Sorption Phenomena", Food and Bioprocess Technology
March 2008 Vol 1, Issue 1, pp 82--90
aw/w = a + b*aw + c*aw^2
a = 1/(Wm*C*K)
b = (C-2)/(Wm*C)
c = K(1-C)/(Wm*C)
parameters
aw: array - water activity
w: array - water content (water mass / total mass) at each water activity point
returns dictionary {Wm: float, C: float, K: float}
"""
y = aw / w
[c, b, a] = np.polyfit(aw, y, 2)
K = (-b + np.sqrt(b**2 - 4*a*c))/(2*a)
C = b / (a*K) + 2
Wm = 1 / (b + 2*K*a)
return {"Wm": Wm, "C": C, "K": K}
In [46]:
params = GAB_params.to_dict("records")[0]
w = GAB(WATER_ACTIVITY, params["Wm"], params["C"], params["K"])
params_calculated = GAB_regress(WATER_ACTIVITY, w)
print ("tabulated values: %s" %params)
print ("regression values: %s" %params_calculated)
In [47]:
def GAB_avg(frac, params):
"""
Calculate the GAB parameters for a composition of components.
parameters:
frac: list of floats - list of mass fractions of the individual components [f_i] ; for i=1...N
params: list of dictionaries - list of GAB parameters [{wm, C, K}_i] ; for i=1...N
returns dictionary {Wm: float, C: float, K: float}
"""
aw = np.arange(0.1, 0.8, 0.05)
w = np.array([mass_water_solid_mixture(aw_i, frac, GAB_params.to_dict("records")) for aw_i in aw])
return GAB_regress(aw, w)
In [48]:
avg = GAB_avg(frac, GAB_params.to_dict("records"))
temp = GAB_params.append(pd.DataFrame(avg, index=["AVERAGE"]))
print (temp)
print (pd.DataFrame(data=frac, index=GAB_params.index, columns=["fraction"]))
water = np.array([GAB(WATER_ACTIVITY, *temp.loc[c].values) for c in temp.index.values])
plt.plot(WATER_ACTIVITY, water.T)
plt.legend(temp.index.values, loc="best")
Out[48]:
In [1]:
def GAB_system(frac, params, m, V, T):
"""
Calculate the GAB parameters for a system of both solids and water vapor.
parameters:
frac: list of floats - list of mass fraction of the individual components [m_i] ; for i=1...N
params: list of dictionaries - list of GAB parameters [{wm, C, K}_i] ; for i=1...N
m: float - total mass of solid (kg)
V: float - Volume of water vapor (m^3)
T: float - Temperature (K)
returns dictionary {Wm: float, C: float, K: float}
"""
aw = np.arange(0.1, 0.8, 0.05)
w = np.array([m * mass_water_solid_mixture(aw_i, frac, GAB_params.to_dict("records"))
+ mass_water_vaport(T, aw_i, V) for aw_i in aw])
return GAB_regress(aw, w)
In [59]:
def water_activity_GAB(w, Wm, C, K):
"""
Calculate the water activity, aw, from the known water content, w, in a substance
and known GAB parameters, Wm, C, K.
From "GAB Generalized Equation for Sorption Phenomena", Food and Bioprocess Technology
March 2008 Vol 1, Issue 1, pp 82--90
aw/w = a + b*aw + c*aw^2 --> c*aw^2 + (b-1/w)*aw + a = 0
solution from quadratic equation
aw = (-(b-1/w) +/- sqrt((b-1/w)^2 - 4*c*a)) / (2*c)
where
a = 1/(Wm*C*K)
b = (C-2)/(Wm*C)
c = K(1-C)/(Wm*C)
parameters
w: float - water content in component (mass water / total mass solid)
Wm: float - GAB parameter
C: float - GAB parameter
K: float - GAB parameter
returns float - water activity
"""
a = 1/(Wm*C*K)
b = (C-2)/(Wm*C)
c = K*(1-C)/(Wm*C)
#arg = np.sqrt((b-1/w)**2 - 4*c*a)
# TODO How do we know which root to use?
#hi = (-1*(b-1/w) + arg) / (2*C)
#lo = (-1*(b-1/w) - arg) / (2*c)
#print lo, hi
# Relationship from Gary"s method (no reference)
aw = (np.sqrt(C) * np.sqrt(C*Wm**2 - 2*C*Wm*w + C*w**2 + 4*Wm*w) - C*Wm + C*w - 2*w) / (2*(C-1)*K*w)
return aw # TODO select the calculation to use
In [60]:
water = GAB(WATER_ACTIVITY, *GAB_params.iloc[0].values)
aw_calculated = np.array([water_activity_GAB(w, *GAB_params.iloc[0].values) for w in water])
plt.plot(WATER_ACTIVITY, aw_calculated)
Out[60]:
In [117]:
def transport_coef_calc(T, m, b, l=1, A=1, f1=1, f2=0, LDPE=False):
"""
Calculate the water vaport transport coefficient of the package at the specified temperature.
parameters
T: float - temperature (K)
m: float - permeability constant
b: float - permeability constant
l: float - thickness of package (m)
A: float - Area of package (m^2)
f1: float - seal constant (unitless)
f2: float - seal constant (mg/day/delta RH)
return float - permeability kg/s/delta RH
"""
T = T - 273.0
if LDPE:
l = l * 39370.1
p = np.exp(b + m/T) / (l/4.)
p = 1.e-3 / (24.*3600) * p
else:
p = np.exp(b + m/T) * f1 + f2
p = 1.e-6 / (24.*3600) * p
return p
In [62]:
def equilibrium(mw, Vg, T, ms, GAB):
"""
Calculate the equilibrium relative humidity (rh) from the partition of water between
a solid and vapor in a closed container (i.e. constant water mass).
mw: float - total mass of water in the system (kg)
Vg: float - volume of the gas phase (m^3)
T: float - temperature (K)
ms: float - mass of solid (kg)
GAB: dictionary - GAB parameters {Wm, C, K}
returns float equilibrium relative humidity
"""
m_sat = mass_density_sat(T) * Vg
def objective(mws):
"""
Objective function for calculating water partitioning. The objective function is equal to the
difference between the relative humidity calculated from the water vapor mass and
the relative humidity calculated from the water in the solid with isotherm provided
by the GAB parameters. Objective is to be minimized by varying the mass of water in
the solid, mws.
Function:
(rh_v - rh_s)**2
where,
rh_v = mwv/m_sat (assuming ideality = p_i/p_i^sat = rh)
mwv = mw - mws
rh_s = water_activity_GAB(mws/ms, Wm, C, K) * ms
"""
return ((mw-mws)/m_sat - water_activity_GAB(mws/ms, GAB["Wm"], GAB["C"], GAB["K"]))**2
res = minimize_scalar(objective)
return res
In [63]:
Vg = 1e-4
T = 323.
ms = 0.1
m_sat = mass_density_sat(T) * Vg
mwv = np.array([mass_water_vapor(T, rh_i, Vg) for rh_i in WATER_ACTIVITY])
mws = ms * np.array([GAB(rh_i, *GAB_params.iloc[0].values) for rh_i in WATER_ACTIVITY])
mw = mwv + mws
mws_calc = np.array([equilibrium(mw_i, Vg, T, ms, GAB_params.to_dict("records")[0]).x for mw_i in mw])
plt.plot(mws, mws_calc)
Out[63]:
In [64]:
def step(dt, rh_o, rh_i, T, Vg_o, Vg_i, p, ms, GAB):
"""
Calculate the incremental moisture transport through a package and subsequent
re-equilbration of water between the vapor phase and the solid phase on the inside of the package.
parameters
dt: float - time increment (s)
rh_o: float - relative humidity outside the package (unitless)
rh_i: float - relative humidity inside the package (unitless)
T: float - temperature (K)
Vg_o: float - volume of headspace in the outer package (m^3)
Vg_i: float - volume of headspace in the inner package (m^3)
ms_o: float - mass of solids in the outer package (kg)
ms_i: float - mass of solids in the inner package (kg)
GAB: dictionary - dictionary of average GAB parameters
returns float - new rh of the inner package
"""
mw = mass_water_vapor(T, rh_i, Vg) + ms * GAB(*p["GAB"])
dmw = p (rh_o - rh_i) * dt
mw_new = mw + dmw
rh_new = equilibrium(mw_new, Vg, T, ms, GAB)
return rh_new
In [16]:
inner_desiccant_GAB = GAB_desiccants.loc[system["inner"]["desiccant"]["type"]][["Wm", "C", "K"]]
outer_desiccant_GAB = GAB_desiccants.loc[system["outer"]["desiccant"]["type"]][["Wm", "C", "K"]]
inner_component_GAB = GAB_components.loc[[component["name"] for component in
system["inner"]["product"]["components"]]][["Wm", "C", "K"]]
inner_component_frac = [component["frac"] for component in system["inner"]["product"]["components"]]
In [ ]:
[step(t, 0.6, 0.1, 323, 1e-3, )]